 function [phi,phi_x,phi_y,phi_xx,phi_xy,phi_yy] = fe_Argyris(x, y)  
% function [phi,phi_x,phi_y,phi_xx,phi_xy,phi_yy] = fe_Argyris(x, y)  
% The evaluation of tge Argyris elements
%
%    x and y is the area coordinates, in different order with xi and eta 
%    2
%    |  \
%    |    \
%    |      \
%    3 ------ 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% for ease of access
x2 = x.*x;
x3 = x2.*x;
x4 = x2.*x2;
x5 = x2.*x3;
y2 = y.*y;
y3 = y2.*y;
y4 = y2.*y2;
y5 = y2.*y3;

phi(:, 1) = 1 - 10*x3 - 10*y3 + 15*x4 - 30*x2.*y2 + 15*y4 - 6*x5 + ...
            30*x3.*y2 + 30*x2.*y3 - 6*y5;
phi(:, 2) = x - 6*x3 - 11*x.*y2 + 8*x4 + 10*x2.*y2 + 18*x.*y3 - ...
            3*x5 + x3.*y2 - 10*x2.*y3 - 8*x.*y4;
phi(:, 3) = y - 11*x2.*y - 6*y3 + 18*x3.*y + 10*x2.*y2 + ...
            8*y4 - 8*x4.*y - 10*x3.*y2 + x2.*y3 - 3*y5;
phi(:, 4) = 0.5*x2 - 1.5*x3 + 1.5*x4 - 1.5*x2.*y2 - 0.5*x5 + 1.5*x3.*y2 + x2.*y3;
phi(:, 5) = x.*y - 4*x2.*y - 4*x.*y2 + 5*x3.*y + 10*x2.*y2 + 5*x.*y3 - 2*x4.*y - 6*x3.*y2 - 6*x2.*y3 - 2*x*y4;
phi(:, 6) = 0.5*y2 - 1.5*y3 - 1.5*x2.*y2 + 1.5*y4 + x3*y2 + 1.5*x2*y3 - 0.5*y5;
phi(:, 7) = 10*x3 - 15*x4 + 15*x2.*y2 + 6*x5 - 15*x3.*y2 - 15*x2.*y3;
phi(:, 8) = -4*x3 + 7*x4 - 3.5*x2.*y2 - 3*x5 + 3.5*x3*y2 + 3.5*x2*y3;
phi(:, 9) = -5*x2.*y + 14*x3.*y + 18.5*x2.*y2 - 8*x4.*y - 18.5*x3.*y2 - 13.5*x2.*y3;
phi(:,10) = 0.5*x3 - x4 + 0.25*x2.*y2 + 0.5*x5 - 0.25*x3.*y2 - 0.25*x2.*y3;
phi(:,11) = x2.*y - 3*x3.*y - 3.5*x2.*y2 + 2*x4.*y + 3.5*x3*y2 + 2.5*x2*y3;
phi(:,12) = 1.25*x2.*y2 - 0.75*x3.*y2 - 1.25*x2.*y3;
phi(:,13) = 10*y3 + 15*x2.*y2 - 15*y4 - 15*x3.*y2 - 15*x2.*y3 + 6*y5;
phi(:,14) = -5*x.*y2 + 18.5*x.*y2 + 14*x.*y3 - 13.5*x3.*y2 - 18.5*x2.*y3 - 8*x.*y4;
phi(:,15) = -4*y3 - 3.5*x2.*y2 + 7*y4 + 3.5*x3.*y2 + 3.5*x2.*y3 - 3*y5;
phi(:,16) = 1.25*x2.*y2 - 1.25*x3.*y2 - 0.75*x2.*y3;
phi(:,17) = x.*y2 - 3.5*x2.*y2 - 3*x.*y3 + 2.5*x3.*y2 + 3.5*x2.*y3 + 2*x.*y4;
phi(:,18) = 0.5*y3 + 0.25*x2.*y2 - y4 - 0.25*x3.*y2 - 0.25*x2.*y3 + 0.5*y5;
phi(:,19) = sqrt(2) * (-8*x2.*y2 + 8*x3.*y2 + 8*x2.*y3);
phi(:,20) = -16*x.*y2 + 32*x2.*y2 + 32*x.*y3 - 16*x3.*y2 - 32*x2.*y3 - 16*x.*y4;
phi(:,21) = -16*x2.*y + 32*x3.*y + 32*x2.*y2 - 16*x4.*y - 32*x3.*y2 - 16*x2.*y3;

phi_x = zeros(size(x,1),21);
phi_y = phi_x;
phi_xx = phi_x;
phi_xy = phi_x;
phi_yy = phi_x;